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Abstract. In the Coulomb blockade regime of a ballistic quantum dot, the distribution of conductance peak spacings is well 
known to be incorrectly predicted by a single-particle picture; instead, matrix element fluctuations of the residual electronic 
interaction need to be taken into account. In the normalized random-wave model, valid in the semiclassical limit where the 
number of electrons in the dot becomes large, we obtain analytic expressions for the fluctuations of two-body and one- 
body matrix elements. However, these fluctuations may be too small to explain low-temperature experimental data. We have 
examined matrix element fluctuations in realistic chaotic geometries, and shown that at energies of experimental interest these 
fluctuations generically exceed by a factor of about 3-4 the predictions of the random wave model. Even larger fluctuations 
occur in geometries with a mixed chaotic-regular phase space. These results may allow for much better agreement between 
the Hartree-Fock picture and experiment. Among other findings, we show that the distribution of interaction matrix elements 
is strongly non-Gaussian in the parameter range of experimental interest, even in the random wave model. We also find that 
the enhanced fluctuations in realistic geometries cannot be computed using a leading-order semiclassical approach, but may 
be understood in terms of short-time dynamics. 
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INTRODUCTION 

There has been much interest in the mesoscopic properties of quantum dots whose single-particle dynamics are chaotic 
01- The generic fluctuation properties of the single-particle spectrum and wave functions in such dots are usually 
described by random matrix theory (RMT) [2]. However, in almost-isolated dots, electron-electron interactions are 
important and must also be taken into account. The simplest model is the constant interaction (CI) model, in which 
the interaction is taken to be the classical charging energy. Charging energy leads to Coulomb blockade peaks in 
the conductance versus gate voltage. Each peak occurs as the gate voltage is tuned to compensate for the Coulomb 
repulsion and an additional electron tunnels into the dot. For a fixed number of electrons, the CI model is essentially 
a single-particle model and RMT can be used to derive the statistical properties of the conductance peak heights yfl. 
While the CI plus RMT model has explained (at least qualitatively) JH 0, Ht] several observed features of the peak 
height fluctuations [g, |8fl, there have been significant discrepancies with experimental data, in particular regarding 
the peak spacing statistics l9l [Iol[ITl[l2tl . Such discrepancies indicate the importance of interactions beyond charging 
energy. 

A more systematic way of treating electron-electron interactions in chaotic ballistic dots is to expand the interaction 
in a small parameter, the inverse of the Thouless conductance gj ~ kL^ \/N, where k is the Fermi wave number, L is 
the typical linear size of the dot (i.e. L — \/V, where V is the area), and N is the number of electrons in the dot. In the 
limit of large Thouless conductance (equivalently, in the semiclassical or many-electron limit), only a few interaction 
terms survive, constituting the interacting part of the universal Hamiltonian loflill . These universal interaction terms 
include, in addition to charging energy, a constant exchange interaction. The inclusion of an exchange interaction 
has explained the statistics of peak heights at low and moderate temperatures as well as the suppression of the peak 
spacing fluctuations jTslflrjll . However, at low temperatures, the peak spacing distribution remains bimodal even when 
the exchange interaction is included, while none of the experimental distributions are bimodal ll9l [Tll[Tol[T2ll . 

For finite Thouless conductance, residual interactions beyond the universal Hamiltonian must be taken into account. 
In a Hartree-Fock-Koopmans [ 17] approach (assuming the Hartree-Fock single-particle wave functions do not change 
as electrons are added to the dot), the peak spacings can be expressed directly in terms of certain (diagonal) two-body 
interaction matrix elements OH- Sufficiently large fluctuations of these interaction matrix elements can explain the 
absence of bimodality in the peak spacing distribution Il8lll9ll . In a diffusive dot, the variance of the matrix elements 
of the screened Coulomb interaction was shown to behave as A 2 /g^ to leading order in 1/gr 12011 . where the single- 



electron mean level spacing A sets the energy scale. However, dots studied in the experiments are usually ballistic. 

An additional contribution to the peak spacing fluctuations originates in surface charge effects 12011 . In a finite size 
system, screening leads to the accumulation of charge on the surface of the dot. The confining one-body potential is 
then modified upon the addition of an electron to the dot. 

In this paper we will summarize some of our recent results on fluctuations of the two-body interaction matrix 
elements and of the surface charge one-body matrix elements in ballistic dots 12111 . We begin by defining the 
matrix elements of interest, and note that their fluctuations can be expressed to leading order in 1/gr in terms of 
spatial correlations within single-electron wave functions. Berry's random wave conjecture 112211 provides the first 
approximation for these correlations, which (in contrast with the situation for diffusive dots) is geometry independent. 
However, the spatial correlator of wave function intensity obtained from the Berry conjecture is not consistent with the 
normalization requirement of the wave functions [23]. We discuss the importance of normalization corrections to the 
random wave correlator, and show how the variances of the two-body and one-body interaction matrix elements may 
be computed in a normalized random wave model. We also find that the distribution of interaction matrix elements 
may be very far from Gaussian, even in a normalized random wave model where the wave functions are very close to 
Gaussian. 

An interesting quantity that we refrain from discussing here is the covariance of interaction matrix elements, relevant 
for understanding spectral scrambling when several electrons are added to the dot [24]. 

We then proceed to study matrix element fluctuations in actual chaotic systems, using a family of modified quarter- 
stadium billiards as an example. We find strongly enhanced fluctuations in comparison with the normalized random 
wave model. Semiclassical corrections due to bounces from the dot's boundaries lead to only a modest increase in 
the fluctuations, and do not correctly predict the scaling with kL in the experimentally relevant range. Insight into the 
underlying mechanism of fluctuation enhancement is obtained by examining a family of quantum maps. An important 
conclusion is that the expansion in 1/gr, while asymptotically correct, can be problematic in quantifying matrix 
element fluctuations in the regime relevant to experiments. Finally, in the last section we study systems beyond the 
chaotic regime, i.e., billiards dominated by marginally-stable bouncing-ball modes as well as billiards with mixed 
dynamics (partly regular and partly chaotic). 



INTERACTION MATRIX ELEMENTS 

The general two-body interaction matrix element for potential v(f ,r) is given by 

v«j8 ;y8~ j J drd? \f a (?) \j/* p (/) v(/ , ?) y y {r) Y 5 (f), ( 1 ) 

where \]f a , y/p, ... are single-electron orbital wave functions. In practice, screening of the residual electron-electron 
interaction causes the range of the interaction to be much smaller than the size of the dot. Modeling the residual 
interaction as a contact interaction, we obtain v a p.yg = AV J v dr\j/^, (?) l/^(r) y/y{r) Ys(r), where the mean single- 
particle level spacing A sets the energy scale. We note that for a contact interaction, exchange terms have precisely the 
same form and need not be considered separately. 

Three distinct situations must be treated. First, the diagonal matrix elements 

v a p = v a p. ap =AV^£/?|y/ a (?)| 2 |y/ i 3(?)| 2 (2) 

(0£ ^ j3) appear when two electrons are found on the same two orbitals a and j3 before and after the interaction. 
Secondly, we have the double-diagonal matrix elements v aa = v a a;aa, where two electrons are annihilated from one 
orbital and created on the same orbital (for a contact interaction, v aa is simply an inverse participation ratio in position 
space). Finally, we may consider the off-diagonal matrix elements v a pyg where the four orbitals are all distinct. We 
will focus mostly on the diagonal matrix elements v a p , but the other two situations are treated similarly. 

The mean V^p, averaged over all pairs of orbitals or over an ensemble, may always be absorbed into the mean field 
part of the Hamiltonian. Thus, we are primarily interested not in the average, but in the fluctuations. To leading order 
in gj ~ kL, the dominant contribution to the variance arises from correlations between the intensities of a single wave 
function at different points 12411 : 

lf^~ = A 2 V 2 J J drdr>C 2 (r,r>) + --- (3) 



where 

C(r,/) = \y,(r)\ 2 \W)\ 2 - \W)\ 2 \W)\ 2 ■ (4) 
Thus, the problem of two-body matrix elements statistics has been reduced to an apparently simpler problem of 
understanding the statistics of individual single-electron wave functions in a chaotic potential. Terms we have omitted 
in Eq. @ involve spatial correlators between | y/ a (?) | 2 and | y/p (r 1 ) | 2 ; such correlators are essential for computing the 

covariance 8v a p8v a y, but are subleading in the calculation of the variance. 



RANDOM WAVE APPROXIMATION 

For a two-dimensional billiard system, the random wave model implies that a typical chaotic wave function may be 
written locally as a random superposition of plane waves at fixed energy h 2 k 2 /2tn. Adopting the usual normalization 
|y(r)| 2 = 1 /V, we obtain an amplitude correlator y/*(r)y/(r') = yJo(k\r — r'|) and intensity correlator 

C™(r,r / ) = ~J? ) (k\r-r'\), (5) 

where j8 = 1 or 2 represents the presence or absence of time reversal symmetry, i.e. the absence or presence of an 
external magnetic field, respectively. 

The intensity correlator C rw (r,r') is valid to leading order in \r — r'|/L, but becomes problematic when applied to to 
all r, r* in the finite area V. Indeed wave function normalization requires the correlator to vanish on average, 

/ drC{r,r>) =0. (6) 
Jv 

However, C™(r, f) in (0 manifestly does not satisfy the condition ||6}. The reason for this failure is that in the random 
wave model, normalization is satisfied only for the ensemble average, i.e., f v dr\\]/(r)\ 2 = 1, but not for each individual 
random superposition of plane waves \l/(r). 

This deficiency can be corrected by introducing the normalized random wave model, in which each "random" wave 
function is normalized in area V, i.e., 

r°™{7) = V{r)/ ( dr\W)\ 2 - (7) 
Jv 

Assuming the deviation from exact normalization in area V is small for the original random waves y/(r), we apply a 
perturbative scheme I2H1 and obtain 

C nOTm (r/)=e(r/) + o(^), (8) 

where 

C(r/) = CT{7,¥) - I / dr a C™{r,r a ) - 1 / d? a C™ + / / dr a dr h Cr7 a ,r b ) ■ (9) 
V Jv V Jv V L JvJv 

The leading-order normalized correlator C(r,r') was previously derived in Ref. i23ll by adding a weak smooth dis- 
order and using the non-linear supersymmetric sigma model. Higher-order contributions may be computed system- 
atically; the full expression for the normalized two-point correlator C norm (r,r') involves all unnormalized n-point 

correlators (| \j/(r\)\ 2 — y)--- (| W^n)\ 2 — y), starting with the unnormalized three-point correlator that gives rise to 
the (9((fcL)~ 3 / 2 ) correction in Eq. ((HJ. In practice the higher-order terms are small; the approximation C already satis- 
fies (0 and can be proven to imply normalization of individual wave functions [21]. However, n-point correlators will 
be important below, when we discuss the shape of matrix element distributions. 



Two-body matrix elements 

Substituting the normalized random wave correlator C into Eq. (O, we obtain 

Tx - a2 3 /2\ 2 \nkL + b g f A 2 \ 



where the leading \nkL/ (kL) 2 term may be obtained already from the unnormalized correlator C™ and depends only 
on the symmetry class, while the shape-dependent coefficient b g may be easily evaluated by an integral involving 
Bessel functions over the area V of interest. 
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FIGURE 1. The variance of the two-body matrix element v a a versus kL in the real (j3 = 1) random wave model: (a) the solid 
curve is the result of exact numerical simulations; (b) the long-dashed line is the result of substituting the normalized random wave 
correlator C(r, f) into Eq. ^3); (c) the short-dashed line is the result of using the unnormalized correlator C W (F, ?'). In the inset, the 
solid line is the difference between the full result (a) and the approximation (b); the dashed line indicates that omitted terms scale 
as jj^j- Here and in all following figures, the single-particle level spacing A, which sets the overall energy scale, is set to unity. 

In Fig. Q] we show the results for a disk geometry (b g = —0.1) in the range 30 < kL < 70, corresponding roughly 
to the parameter range relevant for experiments (i.e., a few hundred electrons in the dot). The discrepancy between 
the analytic result dTOb and exact numerics is indeed 0(1/ (£L) 3 ), resulting in an accuracy of ~ 10% for the analytic 
formula. In practice, the dependence of b g on the geometry is weak, e.g. deforming the disk to an ellipse with an aspect 
ratio of 16 while keeping the area fixed produces a change of only ~ 5 — 6% in the v a p variance for the experimental 
range of kL. Thus, for all practical purposes, shape effects on the v a p variance can be ignored (at least within the 
normalized random wave model). 

Also, at this subleading order, we must in principle take into account the wave number difference 8k = k a — kp. For 
(k a +kp)/2 3> 8k 3> l/L (i.e. for a wave number difference that is quantum mechanically large but classically small), 



this results merely in a modification of the geometry-dependent coefficient in Eq. ( fTOb : 



b g -(\/3)\n8kL. In 



practice, for reasonable separations 8kL, the consequent reduction in the variance is at most 10%, and may be safely 
ignored compared to the much larger dynamical effects to be discussed later. 

The variance of double-diagonal or off-diagonal matrix elements may be computed similarly. To leading order, 
O (\nkL / (kL) 2 ) , all the results arise from integrating the unnormalized correlator (O and differ only by geometry- 
independent combinatorial prefactors. However, normalization-related subtraction gives different results in the three 



situations, and thus the subleading coefficient b g in ( fTOt must be replaced with b' g or b" g in the case of v aa or v 
respectively. Thus the variance ratios converge to universal shape-independent constants in the kL 
convergence is logarithmically slow. For example, in the presence of time-reversal symmetry, 



• °° limit, but the 



y. 



In kL 



2.15 
ln/fcL 



(11) 



to leading order in 1/lnkL. For a disk, this ratio barely reaches 3 in the kL range of experimental interest. This is our 
first indication that beautiful analytic results valid as kL — > °° may not always have relevance to experiments, even 
when the dot contains hundreds or thousands of electrons. 



One-body matrix elements 

When an electron is added to the finite dot, charge accumulates on the surface and its effect can be described by a 
one-body potential energy ~f(r). The diagonal matrix elements of "f(r) are given by 



y aa = ( dr\xi/ a (r)\ 2 f(r). 
Jv 



(12) 



Again, we wish to express the variance in terms of the wave function intensity correlator: 



Svk = I ( dfdf -r{r)C{r,f)f{f), (13) 
Jv Jv 

Because only one power of C appears in the variance, the integral is dominated by distant pairs of points \r — f] ~ L, 
and scales as 1 / kL: 

^=f£+°(w) • <14) 

where c ? is a shape-dependent dimensionless coefficient. We note that normalization-related subtraction of the corre- 
lator, which had only a moderate effect on the two-body matrix element fluctuations, here reduces the variance by a 
full order of magnitude, resulting in a very small dimensionless coefficient: c g = 0.035 for a disk, and even smaller for 
less symmetric shapes. This is due the fact that the integrand (fT3l l ceases to be everywhere positive after subtraction 
(0, in contrast with the integrand in (0. The analytic expression < TT~4T > is in excellent agreement with exact Monte 
Carlo simulations (not shown). Due to the small size of the coefficient c g , the numerical value of is smaller in the 
physically interesting kL regime than the corresponding result for the two-body matrix element variance 5v^, despite 
the fact that the former is parametrically larger in a 1 /kL expansion. 



Matrix element distributions 

The central limit theorem implies that all interaction matrix elements in the random wave model must be distributed 
as Gaussian random variables as kL — > °°. This justifies our focus so far on the variance (and covariance) of these 
matrix elements, to the neglect of higher moments. However, we have seen above that non-universal finite kL effects 
are sometimes significant in the experimentally relevant regime kL < 70, e.g., for the variance ratio ( fTTT ). Thus, we 
should look explicitly at the shape of matrix element distributions for finite kL. 
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FIGURE 2. The distribution of diagonal interaction matrix elements v a p is shown for real random waves in a disk (dashed 
curve) and for actual eigenstates in a modified quarter-stadium billiard geometry (solid curve, see below) at kL = 70. A Gaussian 
distribution with the same mean and variance as the random wave distribution is shown as a dotted curve for comparison. 

Data for the diagonal matrix elements v a p in shown in Fig. [2] Results for other shapes and for other matrix elements 
(e.g., v a a, v apy8> or v a) are qualitatively similar. We note that the numerically obtained interaction matrix element 
distribution for random waves (dashed curve) has a long tail on the right side as compared with a Gaussian distribution 
of the same mean and variance (dotted curve). In other words, there is an excess of anomalously large matrix elements, 
compensated for by a reduction in the median to a value slightly below A. 

Deviations from a Gaussian shape can be quantified by considering higher moments. For large kL, we may estimate 
these higher moments in a manner analogous to our estimate for the variance in Eq. (O, but using n— point correlators. 
After some calculation i2ll) . we obtain the skewness 




and excess kurtosis 

72 = I; = b *g c 4 

K» } V 

Because the decay is only logarithmic, j\ and 72 always remain > 1 for values of kL relevant in the experiments. The 
same holds true for other matrix elements and for the higher moments. Therefore, matrix element distributions are 
predicted to be strongly non-Gaussian, even within the random wave model. 




CHAOTIC BILLIARDS 

We now investigate how dynamical effects may modify the fluctuations of interaction matrix elements beyond the 
(normalized) random wave model. As an example, we use a modified quarter-stadium billiard geometry B25I1 . where 
the quarter-circle has radius R and the straight edge of length a has been replaced by a parabolic bump to eliminate 
bouncing-ball modes (Fig. [3). The system has been verified numerically to be fully chaotic for the range of parameters 
used. Variation of the bump size s allows us to check the sensitivity of the results to details of the billiard geometry 
while maintaining the chaotic character of the classical dynamics. Furthermore, variation of the parameter a allows us 
to control the degree of classical chaos. 




FIGURE 3. A modified quarter-stadium geometry with parameters a and s is used to illustrate dynamical effects on matrix 
element fluctuations. In the Figure, we set the quarter-circle radius R = 1. The random wave contribution to the correlator, C™{7\ , fj) 
is schematically indicated by a dashed line, and a typical dynamical contribution by a dotted line. 

We first study the variance of the diagonal interaction matrix elements v a p . Typical numerical results are shown 
in Fig. |U We note the large enhancement of the billiard results over the random wave model. To understand this 
discrepancy, we first compare the exact numerical result for 8v 2 a p with the first term on the right hand side of 

Eq. ([3J, where the intensity correlator for normalized random waves C(r,r') is replaced by the billiard correlator 
C blll (r, r') (calculated numerically for the appropriate billiard system). The discrepancy is immediately reduced to a 
~ 5 — 10% level, which is comparable to the 0(1 /(kL) 3 ) discrepancy observed in the random wave model (see inset 
of Fig. Q]). Thus, the large enhancement of v a p fluctuations over the random wave prediction can be traced directly to 
an enhancement in the single-particle intensity correlator C blll (r,r') over C(r,r'). 

Can this dynamical enhancement of the intensity correlator (as compared with a random wave model) be obtained 
using a semiclassical approach? The random wave correlator C™(r,r l ) may be interpreted semiclassically as arising 
from straight-line free propagation, as indicated by the dashed line in Fig. [3] As discussed by several authors 11261 12711 . 
additional contributions to the correlator in a specific dynamical system can be associated with classical trajectories 
that bounce off the boundary n times on their way from r to f , such as the one indicated by a dotted line in the same 
figure. We obtain 12 ill 

*5 ^ v 1J/^^^ 1+0 {w?) <17) 

2 3 (2\ 2 (lnkL + b g )+b sc 
n\p) (kLf 
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FIGURE 4. The enhancement of the v a p variance over the random wave prediction is shown for sevarl billiards with Neumann 
boundary conditions. Solid lines: a = 0.25, s = 0.1, 0.2. Dashed lines: a = 1.0, s = 0.1, 0.2. 

where b sc is a classical constant that in practice must be determined numerically by performing the integral in Eq. STfi . 
Random wave and semiclassical contributions to C blll (r, f) are of the same order except for \r — r'\ <C L; it is these 
short-distance pairs that result in a logarithmic enhancement of the random-wave term. 

The coefficient b sc may in practice be quite large, even for generic chaotic systems, such as the modified stadium 
billiard. In a diagonal approximation, b sc ~ r c j as where r c i as is the decay time of classical correlations. Qualitatively, 
this is consistent with Fig. |4] as fluctuations are observed to be consistently larger for the less chaotic a = 0.25 billiard, 
as compared with the a = 1.00 billiard. We note that both billiards are "generic", in the sense that they are not fine- 
tuned to obtain an anomalously long time scale r c i as . We also note that varying the bump size s has a very weak effect 
on the matrix element statistics (as long as s is large enough to destroy the bouncing-ball modes) and serves instead to 
provide an estimate of the statistical fluctuations. 

However, a closer look at the data shows that the numerical results cannot be explained fully by semiclassical 
arguments, since ( fTST l converges to the random wave limit ( TTOb as kL — > oo, while the enhancement factor in Fig. |4] 
grows with kL in the parameter range of experimental interest. 
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FIGURE 5. The two-body matrix element variance 5 for a quantum map, as a function of Hilbert space dimension N. From top 
to bottom, the three solid lines represent data for dominant orbit stability exponent Ao = 0.25, 0.50, 1.00. The three dashed lines 
indicate the asymptotic 1 /N 2 behavior for each case in the semiclassical regime of large N. 

This anomalous behavior results from a combined effect of two factors: the large numerical value of b sc for generic 
dynamical systems, and saturation of the 1 / (kL) 2 behavior at moderate (< 100) values of kL. As the classical system 
becomes less unstable and the correlation time T co „ increases, b sc also increases, leading to greatly enhanced matrix 
element variance at very large values of kL ( fl~8T >. Because the variance is bounded above independent of kL, the 
(kL)~ 2 growth in the variance necessarily breaks down at moderate values of kL. This saturation sets in at ever larger 
values of kL as the system becomes less unstable. Alternatively, one may note that the natural expansion parameter for 
interaction matrix element fluctuations in a dynamical system is not 1/kL but rather the inverse Thouless conductance 
gj 1 ~ T corr /kL, and the semiclassical contribution with prefactor b sc in Eq. ( fT~8l > is the leading 0(gj 2 ) effect in such an 
expansion. Terms of third and higher order in g^ 1 , although formally subleading and not included in a semiclassical 
calculation, become quantitatively as large as the leading 0(gj 2 ) term when gj falls below some characteristic value. 



Furthermore, if one considers chaotic billiards with a long correlation decay time T corr , the importance of formally 
subleading terms in the gj 1 expansion will extend to quite large values of kL. 

The above assertions are explicitly confirmed for a quantum map model, which has scaling behavior analogous to 
that of a two-dimensional billiard, with the number of states N = 2n/h playing the role of semiclassical parameter 
kL = pL/hin the billiard Il30ll3lll . As in the billiard, a free parameter Ao in the definition of the map allows for control 
of the classical correlation decay time T corr . We see in Fig.[5]that the expected N~ 2 behavior of the variance is observed 
at sufficiently large N, for all three families of maps considered. Furthermore, the prefactor multiplying N~ 2 in each 
case agrees with that obtained from a semiclassical calculation, and as expected this prefactor grows with increasing 
classical correlation time T con - (corresponding to a decrease in the chaoticity of the system). We also see in Fig.[5]that 
even for a "typical" chaotic system (i.e., r corr ~ 1), strong deviations from the 1/N 2 law appear already below N » 80. 
Such deviations extend to even larger N for chaotic systems with slower classical correlation decay. This suggests 
that the large-Af or large-fcL expansion, though theoretically appealing and asymptotically correct, is problematic in 
describing the quantitative behavior of interaction matrix element fluctuations for real chaotic systems in the physically 
interesting energy range. 

The behavior of v aa , v a py§, and v a in real dynamical systems may be studied similarly, again using the normalized 
random wave predictions as the baseline for comparison. The enhancement of the variance at large kL is particularly 
dramatic in the case of double-diagonal matrix element fluctuations. This is consistent with the reasonable expectation 
that dynamical effects lead to particularly strong deviations from random wave behavior in a modest fraction of the 
total set of single-particle states, such as those associated with scarring on unstable periodic orbits |28]. Such deviations 
lead to a significant tail in the v aa distribution, but have a minimal effect on the distribution of off-diagonal matrix 
elements, since it is unlikely for all four wave functions y/ a , y/p, \]/y, and y/§ to be strongly scarred or antiscarred on 
the same orbit. 

We can also go beyond the variance to investigate higher moments of the matrix element distribution for actual 
chaotic systems. A typical distribution for diagonal two-body matrix elements v a p in a modified quarter-stadium 
billiard with a = 0.25 and s = 0.1 is shown in Fig.|2](solid line). Since the approach to Gaussian behavior is already 
very slow in the case of random waves, it is not surprising to find even stronger deviations from a Gaussian shape for 
matrix elements in real chaotic systems at the same energies. For example, the excess kurtosis ji increases from 8.3 at 
kL = 70 to 20.9 at kL = 140, while dropping from 3.7 to 3.3 in the random wave model. Similar behavior is obtained 
for the skewness and for other matrix elements. Clearly, the distribution tails are very long, and the assumption of 
Gaussian matrix element distributions is even less justified for real chaotic systems than it was in the random wave 
model. 



BEYOND THE CHAOTIC REGIME 

Finally, we consider fluctuations of matrix elements in systems that are not fully chaotic. Here no universal behavior 
is expected but we shall see that in such systems the variance can be enhanced much more than in fully chaotic 
systems ll29ll . We use the modified quarter-stadium billiard with s = or a < 0. The choice s = corresponds to the 
original Bunimovich stadium, whose quantum fluctuation properties are dominated by the marginally-stable bouncing- 
ball modes, while a < corresponds to a lemon billiard, which has a classically mixed, or soft chaotic, phase space. 

In contrast with the In kL/(kL) 2 falloff in the v a p variance predicted for fully chaotic dynamics by Eq. (fT8l. i n 
the case of regular or mixed dynamics we expect fcL-independent matrix element fluctuations of order unity 12111 . 
Equivalently, we have very large enhancement, scaling as (kL) 2 /InfcL, of the matrix element variance in mixed 
dynamical systems, over the random wave prediction. 

The diagonal matrix element variance 5v 2 ^ is computed as a function of kL for two typical mixed phase-space 
quarter-lemon billiards and compared with the normalized random wave prediction (dashed lines in Fig. |6). As 
expected, the enhancement becomes more pronounced at larger kL. Enhancement of an order of magnitude or 
more over random wave behavior can easily be obtained for physically interesting values of kL. The most dramatic 
enhancement is observed for the a — —0.25 quarter-lemon billiard, which is closer to integrability. 

Behavior intermediate between hard chaos and mixed chaotic/regular phase space is obtained in the presence of 
families of marginally stable classical trajectories, such as the "bouncing ball" orbits of the stadium billiard. In 
the quarter stadium billiard (s = in Fig. [3]), exceptional states associated with such orbits are concentrated in the 
rectangular region of the billiard and constitute a fraction ~ 1 /(kL) x l 2 of the total set of states 1331 . When a and j3 are 
both bouncing ball states, 8v a p ~ v^Tj" ~ A, just as would be the case for regular states concentrated in a finite fraction 
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FIGURE 6. Enhancement of v a p variance as compared with the random wave prediction for a = 0.25, 1.00 quarter-stadium 
billiards (solid lines); a = —0.25, —0.50 quarter-lemon billiards (dashed lines). 



of the available coordinate space. These special matrix elements dominate the variance, leading to 



Numerical data for quarter-stadium billiards is shown by solid lines in Fig. [6] The stronger fluctuations are observed 
in the less chaotic a = 0.25 stadium. 

Again, analogous results are obtained for the other matrix elements (not shown). 



SUMMARY 

We have studied fluctuations of two-body and one-body matrix elements in ballistic quantum dots as a function of 
semiclassical parameter kL. Understanding the quantitative behavior of these fluctuations is important for the proper 
analysis of peak spacing statistics in the Coulomb blockade regime. 

We find that the variance and higher cumulants of two-body and one-body matrix elements may be expressed 
in terms of spatial correlations within single-particle Hartree-Fock wave functions. For the purpose of computing 
the variance of two-body matrix elements v aa , v a p, and v a pyg, these correlations may be approximated, to leading 
order in kL for a chaotic system, by correlations given by Berry's random wave model. The calculation results in a 
variance scaling as \nkL/(kL) 2 , with universal prefactors depending only on the symmetry class of the system. Shape- 
dependent effects on the variance enter at 0(1/ (kL) 2 ), where the random wave intensity correlator must be corrected 
to satisfy individual wave function normalization in finite volume. In the normalized random wave model, ratios such 
as ^ v aa/ ^ y a/3 conver g e oni y with a logarithmic rate in the kL^ °° limit; as a result, the asymptotic values of such 
ratios are of little relevance in the regime of experimental interest. 

The variance of one-body matrix elements v a is affected by normalization even at leading order, resulting in 
0(1 /kL) scaling in a random wave model, with a shape-dependent prefactor. Both two-body and one-body matrix 
elements follow, already within the random wave model, a strongly non-Gaussian distribution, for all physically rea- 
sonable values of kL. Thus, higher cumulants of these matrix elements will be important in peak spacing calculations, 
especially in the case of two-body matrix elements where we show that the approach to a Gaussian distribution is 
logarithmically slow. 

Dynamical effects, associated with nonrandom short-time behavior in actual chaotic systems, are formally sublead- 
ing for two-body matrix elements, and of the same order as the random wave prediction for one-body matrix elements. 
In practice, however, these effects can easily lead to enhancement by a factor of 3 or 4 of the variance in both one-body 
and two-body matrix elements, for experimentally relevant values of kL in reasonable hard chaotic geometries. The 
size of these effects scales in each case as a power of r c i as , a time scale associated with approach to ergodicity in the 
corresponding classical dynamics. Random wave behavior is recovered as r c i as — > 0. In typical geometries, dynamical 
effects on matrix element fluctuations cannot be properly computed in a semiclassical approximation, as higher-order 
terms are quantitatively of the same size as the leading-order semiclassical expression in the kL range of experimen- 
tal interest. The approach to semiclassical scaling at very large values of kL as well as saturation of matrix element 
fluctuations at moderate to small kL are investigated in the context of a quantum map model. 



Systems with a mixed chaotic-regular phase space or with families of marginally stable classical orbits show 
even stronger enhancement of matrix element fluctuations as compared with the random wave model. The expected 
asymptotic scaling with kL of the fluctuations in these cases is discussed, and is very different from the scaling found 
in chaotic systems. 

Our calculations strongly indicate that statistics of actual chaotic single-particle systems, including dynamical 
effects, are needed to make a proper comparison between Hartree-Fock-Koopmans theory and experiment. A better 
understanding of correlations in single-particle systems is then essential to compute observable properties of the 
interacting many-electron system. Furthermore, these correlations need to be understood beyond the naive leading 
order semiclassical approximation, to allow comparison with experiments, which are generally performed at moderate 
values of the semiclassical parameter. 
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